################################################################################
########################     Analysis and Plots     ############################
################################################################################

library(Matching)
library(tidyverse)

################################################################################

rm(list=ls())
load("~/data/virginclip_pop.RData")

################################################################################

# 1. Dryland: peat depth = 0

virginclip_pop <- virginclip_pop %>%
  filter(peat_depth_GW==0)

# did outcome vars

virginclip_pop$did2016 <- virginclip_pop$forestha_2016 - virginclip_pop$forestha_2011
virginclip_pop$did2017 <- virginclip_pop$forestha_2017 - virginclip_pop$forestha_2011
virginclip_pop$did2018 <- virginclip_pop$forestha_2018 - virginclip_pop$forestha_2011

# time indices

t1 = 2010
t11 = 2011
t0 = 2004

t2 = 2016
t3 = 2017
t4 = 2018

# ddd outcome vars

virginclip_pop$didpre = virginclip_pop$forestha_2010 - virginclip_pop$forestha_2004

virginclip_pop$ddd2016 <- virginclip_pop$did2016 - ((t2-t11)/(t1-t0))*virginclip_pop$didpre
virginclip_pop$ddd2017 <- virginclip_pop$did2017 - ((t3-t11)/(t1-t0))*virginclip_pop$didpre
virginclip_pop$ddd2018 <- virginclip_pop$did2018 - ((t4-t11)/(t1-t0))*virginclip_pop$didpre

# ps model
ps_virginclip_pop <- glm(moratorium ~ forestha_2005 + forestha_2006 + forestha_2007 + 
                       forestha_2008 + forestha_2009 + forestha_2010 +  elev + slope + 
                       distroads + dist_cities + dist_palms_km + dist_timber_km + dist_logging_km +
                       agc_bcc_mean + pop_2005 + pop_2010, family = binomial, data = virginclip_pop)
summary(ps_virginclip_pop)

# attach propensity score to data
virginclip_pop$pscore<- ps_virginclip_pop$fitted

###### DID REGRESSIONS ######

PSM_DID_virginclip_pop2016 <- Match( Y = virginclip_pop$did2016, Tr = virginclip_pop$moratorium,
                                 X = virginclip_pop$pscore, estimand = "ATT", M = 1, ties = T,
                                 replace = T, distance.tolerance = 1e-15, caliper = 0.0001)

summary(PSM_DID_virginclip_pop2016)


PSM_DID_virginclip_pop2017 <- Match( Y = virginclip_pop$did2017, Tr = virginclip_pop$moratorium,
                                 X = virginclip_pop$pscore, estimand = "ATT", M = 1, ties = T,
                                 replace = T, distance.tolerance = 1e-15, caliper = 0.0001)

summary(PSM_DID_virginclip_pop2017)

PSM_DID_virginclip_pop2018 <- Match( Y = virginclip_pop$did2018, Tr = virginclip_pop$moratorium,
                                 X = virginclip_pop$pscore, estimand = "ATT", M = 1, ties = T,
                                 replace = T, distance.tolerance = 1e-15, caliper = 0.0001)

summary(PSM_DID_virginclip_pop2018)



###### DDD REGRESSIONS ######


PSM_DDD_virginclip_pop2016 <- Match( Y = virginclip_pop$ddd2016, Tr = virginclip_pop$moratorium,
                                 X = virginclip_pop$pscore, estimand = "ATT", M = 1, ties = T,
                                 replace = T, distance.tolerance = 1e-15, caliper = 0.0001)

summary(PSM_DDD_virginclip_pop2016)

PSM_DDD_virginclip_pop2017 <- Match( Y = virginclip_pop$ddd2017, Tr = virginclip_pop$moratorium,
                                 X = virginclip_pop$pscore, estimand = "ATT", M = 1, ties = T,
                                 replace = T, distance.tolerance = 1e-15, caliper = 0.0001)

summary(PSM_DDD_virginclip_pop2017)


PSM_DDD_virginclip_pop2018 <- Match( Y = virginclip_pop$ddd2018, Tr = virginclip_pop$moratorium,
                                 X = virginclip_pop$pscore, estimand = "ATT", M = 1, ties = T,
                                 replace = T, distance.tolerance = 1e-15, caliper = 0.0001)

summary(PSM_DDD_virginclip_pop2018)

# remove data and save results
rm(virginclip_pop, ps_virginclip_pop)

save.image("~/results/dryland-np-did-ddd-161718.RData")


### Peatland ###

rm(list=ls())
load("/users/sileci/np-ddd//virginclip_pop.RData")

################################################################################

# 2. Peatland: peat depth > 0

virginclip_pop <- virginclip_pop %>%
  filter(peat_depth_GW>0)

# did outcome vars

virginclip_pop$did2016 <- virginclip_pop$forestha_2016 - virginclip_pop$forestha_2011
virginclip_pop$did2017 <- virginclip_pop$forestha_2017 - virginclip_pop$forestha_2011


# time indices

t1 = 2010
t11 = 2011
t0 = 2004

t2 = 2016
t3 = 2017

# ddd outcome vars
virginclip_pop$didpre = virginclip_pop$forestha_2010 - virginclip_pop$forestha_2004

virginclip_pop$ddd2016 <- virginclip_pop$did2016 - ((t2-t11)/(t1-t0))*virginclip_pop$didpre
virginclip_pop$ddd2017 <- virginclip_pop$did2017 - ((t3-t11)/(t1-t0))*virginclip_pop$didpre


# ps model
ps_virginclip_pop <- glm(moratorium ~ forestha_2005 + forestha_2006 + forestha_2007 + 
                       forestha_2008 + forestha_2009 + forestha_2010 +  elev + slope + 
                       distroads + dist_cities + dist_palms_km + dist_timber_km + dist_logging_km +
                       agc_bcc_mean + pop_2005 + pop_2010, family = binomial, data = virginclip_pop)
summary(ps_virginclip_pop)

# attach propensity score to data
virginclip_pop$pscore<- ps_virginclip_pop$fitted

###### DID REGRESSIONS ######

PSM_DID_virginclip_pop2016 <- Match( Y = virginclip_pop$did2016, Tr = virginclip_pop$moratorium,
                                 X = virginclip_pop$pscore, estimand = "ATT", M = 1, ties = T,
                                 replace = T, distance.tolerance = 1e-15, caliper = 0.0001)

summary(PSM_DID_virginclip_pop2016)



PSM_DID_virginclip_pop2017 <- Match( Y = virginclip_pop$did2017, Tr = virginclip_pop$moratorium,
                                 X = virginclip_pop$pscore, estimand = "ATT", M = 1, ties = T,
                                 replace = T, distance.tolerance = 1e-15, caliper = 0.0001)

summary(PSM_DID_virginclip_pop2017)


###### DDD REGRESSIONS ######

PSM_DDD_virginclip_pop2016 <- Match( Y = virginclip_pop$ddd2016, Tr = virginclip_pop$moratorium,
                                 X = virginclip_pop$pscore, estimand = "ATT", M = 1, ties = T,
                                 replace = T, distance.tolerance = 1e-15, caliper = 0.0001)

summary(PSM_DDD_virginclip_pop2016)

PSM_DDD_virginclip_pop2017 <- Match( Y = virginclip_pop$ddd2017, Tr = virginclip_pop$moratorium,
                                 X = virginclip_pop$pscore, estimand = "ATT", M = 1, ties = T,
                                 replace = T, distance.tolerance = 1e-15, caliper = 0.0001)

summary(PSM_DDD_virginclip_pop2017)


# remove data and save results

rm(virginclip_pop, ps_virginclip_pop)
save.image("~/results/peatland-np-did-ddd-1617.RData")

